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Controlling the onset of traveling pulses in excitable media 
by nonlocal spatial coupling and time-delayed feedback 
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The onset of pulse propagation is studied in a reaction-difltusion (RD) model with control by 
augmented transmission capability that is provided either along nonlocal spatial coupling or by time- 
delayed feedback. We show that traveling pulses occur primarily as solutions to the RD equations 
while augmented transmission changes excitability. For certain ranges of the parameter settings, 
defined as weak susceptibility and moderate control, respectively, the hybrid model can be mapped to 
the original RD model. This results in an efi'ective change of RD parameters controlled by augmented 
transmission. Outside moderate control parameter settings new patterns are obtained, for example 
step-wise propagation due to delay-induced oscillations. Augmented transmission constitutes a 
signaling system complementary to the classical RD mechanism of pattern formation. Our hybrid 
model combines the two major signalling systems in the brain, namely volume transmission and 
synaptic transmission. Our results provide insights into the spread and control of pathological 
pulses in the brain. 



Traveling pulses are of fundamental importance 
in neuroscience. They propagate not only infor- 
mation along nerve fibers, but are also related to 
pathological phenomena. Examples are cell depo- 
larizations that lead to a temporary complete loss 
of normal cell functions in migraine and stroke. 
This state spreads in cortical tissue via chemi- 
cal signals that diffuse through the extracellular 
space. We study these spreading depolarization 
pulses in a standard react ion- diffusion model and 
suggest to augment the transmission capabilities 
such that they reflect the cortical structural and 
functional connectivity. With this new modal- 
ity, we investigate control of emerging spread of 
pathological states in the brain. 



1. INTRODUCTION 

Within the last years control of complex dynamics has 
evolved as one of the central issues in applied nonlin- 
ear science Major progress has been made in neu- 
roscience, among other areas, by extending methods of 
chaos control, in particular time-delayed feedback Q , to 
spatio-temporal patterns d, |3j [H| , and by developing ap- 
plications in the field of biomedical engineering [y, |7|. 
In this study, control is introduced to suppress spatio- 
temporal pattern formation. Our emphasis is on under- 
standing the recruitment of cortical tissue into dysfunc- 
tional states by traveling pulses of pathological activity, 
in particular, on internal cortical circuits that provide 
augmented transmission capabilities and that can pre- 
vent such events. Our long-term aim is to design strate- 
gies that either support the internal cortical control or 
mimic its behavior by external control loops and trans- 
late these methods into applications. 

There is growing experimental evidence that particu- 
lar spatio-temporal pulse patterns in the human cortex. 



called cortical spreading depression, cause transient neu- 
rological symptoms during migraine [1, Q . Similar pulses 
occur after stroke, called periinfarct depolarization, and 
contribute to the loss of potentially salvageable tissue, 
i. e., tissue at risk of infarction These dysfunctional 
states of the cortex are also referred to as spreading depo- 
larizations (SD) to point out the nearly complete depo- 
larizations of cortical cells and its spread as the common 
aspect in these patterns. 

SD is usually called a cortical wave, not pulse, which 
might cause some confusion. In many cases these two 
terms can be used interchangeably. We adhere to a pre- 
cise mathematical terminology, referring to a traveling 
pulse as a localized event with a spatial profile having 
a single front and back while a wave usually refers to a 
periodic spatial profile. SD is a pulse in this terminology, 
the strict use of which is necessary because the spatio- 
temporal pattern formation of SD occurs in cortical tissue 
being only weakly susceptible to SD [ll|. The definition 
of weak susceptibility (Sec. [2]) depends on a bifurcation 
for which one must strictly distinguish between solitary 
and periodic wave forms. 

The pulse of SD extends in the cortex over several cen- 
timeters with a remarkably slow speed of several mil- 
limeter per minute. Accordingly, the mathematical de- 
scription of SD considers large-scale neuronal activity in 
populations of neurons rather than ion channels in the 
membrane of nerve fibers, and the spatial coupling is 
provided by volume transmission, which is essentially a 
diffusion process. The first mathematical model of SD 
was proposed by Hodgkin and Grafstein [l^ based on 
bistable dynamics and diffusion. This model, 
however, describes only the front dynamics of SD. We 
suggest (Sec. [3]) two extensions to this model to study 
the onset of SD. The first is a necessary extension to ob- 
tain onset behavior in the Hodgkin-Grafstein model. We 
include in a generic form a recovery process for the pulse 
trailing edge. In a second step, we extend this model fur- 
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ther by nonlocal and time-delayed signal transmission to 
study the effect of internal control provided as a feedback 
to the traveling pulse. Results and conclusions are given 
in Sec. m and Sec. El respectively. 



2. CLASSIFICATIONS OF EXCITABILITY IN 
LOCAL ELEMENTS AND EXTENDED MEDIA 

Excitability, as a property of a single element, is based 
on threshold behavior and therefore requires a nonlin- 
ear process with a stable fixed point. If the system is 
sufficiently perturbed from this fixed point, it returns af- 
ter a large excursion in phase space, emitting a spike 
[l^. If excitable elements are locally interconnected, a 
new behavior can emerge, namely the capacity to prop- 
agate a sustained pulse through this spatially extended 
system. This emergent property defines a medium as be- 
ing excitable, also termed an active medium. Excitable 
elements and excitable media differ in their response to 
superthreshold stimulation. The response of an excitable 
element to a superthreshold stimulation will eventually 
end in the stable fixed point value of the steady state 
of this clement. In contrast, the response of an ex- 
citable medium, which is initially in the homogeneous 
steady state, to a superthreshold stimulation results in 
approaching a new attractor, the pulse solution. This 
different behavior of local elements and spatial media in- 
dicates that there are also different ways to classify ex- 
citability in elements and media. 

In Sec. 12.11 definitions arc provided for excitability of 
local elements with a focus on neural systems. Some of 
our results concerning control of spatial-temporal pat- 
terns can be explained by considering the effect of con- 
trol on the local dynamics fSec. I4.2~T|) . Furthermore, we 
will consider local dynamics with the aim to extend the 
Hodgkin-Grafstein model (Sec. [3]). However, our main fo- 
cus is on spatial excitability (Sec. 12. 2p and its control by 
nonlocal and time-delayed feedback described in Sec. 14.11 
and Sec. 14.21 respectively. For a thorough treatment of 
local dynamics, in particular in neural systems, see for ex- 
ample Ref. [31 or for more complex discharge patterns, 
like bursting, Ref. p^ . 

2.1. Local excitability 

A generic mechanism of local excitability requires a 
certain configuration of trajectories in the phase space 
of a single excitable element (inset of Fig. [1]). This con- 
figuration usually results from the parameter vicinity of 
an oscillatory regime whose large amplitude limit cycle 
is suddenly destructed [ll] . 

The parameter space is schematically depicted in 
Fig. [TJ where the white area corresponds to the oscilla- 
tory regime, and the colored areas mark various regimes 
of excitability and non-excitability. If parameters set- 
tings are in the excitable regime (to the right of the thick 



blue dashed line), a rest state (fixed point) is the only 
attractor. There exists also a threshold close to this rest 
state, e.g. a "sharp" separatrix or trajectory in phase 
space (thin dashed line, in the inset). Trajectories start- 
ing on the near side of the threshold (green dotted) ap- 
proach the rest state directly (subthreshold response), 
while those starting on the far side (red solid) perform a 
large excursion in phase space, guided by the ghost of a 
destructed large amplitude limit cycle, before returning 
to the rest state (superthreshold response). Note that we 
refer to excitable elements and not to neurons because 
also a population of neurons described in a firing rate 
model can behave like an excitable element. 

The way in which, by changing a bifurcation parame- 
ter, the dynamics of the local element changes from ex- 
citable to oscillatory behavior, that is, the way the large 
amplitude limit cycle is built up, is used to classify the 
type of excitability in a single excitable element. This can 




FIG. 1: Parameter space of excitable systems illustrating 
a universal scheme for classifications of both local and spa- 
tial excitability. Bifurcation lines for excitable elements of an 
active media are marked by thick blue (dashed) lines. Tran- 
sitions from the excitable (magenta, red, orange) to the oscil- 
latory regime (white) occur usually either through a saddle- 
node infinite period bifurcation (SNIPER) causing excitabil- 
ity of type I, or via a Hopf bifurcation (H) causing excitability 
of type 11. Between the SNIPER and a further saddle-node 
(SN) bifurcation line (magenta) three fixed points (stable 
node, unstable focus, and saddle) exist in the local excitable 
elements. The thick white solid lines dC, dM, and dR mark 
bifurcations of active media where the spatio-temporal pat- 
tern in 2D changes (complex, meandering, and rigid rotating 
spiral patterns, to the left of dC, dM and dR, respectively). 
At dP the medium becomes non-excitable. Sustained pulses 
do not exist in the yellow regime but there is a transient 
propagating pulse. To the right of dS (green) the transient 
activation radius becomes zero. The insets (thin lines) indi- 
cate the corresponding configuration of trajectories in phase 
space, with limit cycle (LC), sharp separatrix (S), and fixed 
point (FP). 
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FIG. 2: (Color online) Scheme of the classification of ex- 
citability according to spatio-temporal patterns (spiral-wave, 
particle-like, retracting, collapsing, no spread) in 2D: wave 
front at two instances in time {ti: black, t2: gray) and tra- 
jectory of open wave ends (dashed). This scheme illustrates 
the effect of control on the propagation boundary (dP) by ad- 
ditional nonlocal and time-delayed transmission capabilities. 
dP is controlled by shifting the bifurcation in parameter space 
(horizontal (multi-colored) line). This affects the regime of 
weak susceptibility centered around dP and bounded by the 
rotor boundary OR and the spreading boundary dS. 



happen in different ways, two of which are usually distin- 
guished [3] . Type I excitability obtains its characteristic 
features from a saddle-node infinite period (SNIPER) bi- 
furcation: the frequency of the emerging periodic orbit 
tends to zero, while the amplitude starts with a finite 
value. Type II excitability is caused by a Hopf bifur- 
cation, which is characterized by the features that the 
periodic orbit emerges with zero amplitude and nonzero 
frequency. In practice, the main characteristic for type II 
excitability is, however, the onset frequency, because for 
this type a canard explosion renders the zero amplitude 
practically invisible. 



2.2. Onset of pulse propagation in spatial 
excitability 

The important criterion for spatial excitability is the 
existence of a traveling pulse solution. In contrast to an 
excitable local element and its corresponding phase space 
configuration (inset in Fig. [T]), excitability in active me- 
dia is based on bistability. A superthreshold stimulation 
takes the system from the homogeneous steady state into 
the basin of attraction of the pulse solution. The classi- 
fication of spatial excitability in active media is based on 
this configuration in phase space. For example, the pri- 
mary rough classification into excitable and non-excitable 
media is based on the existence and non-existence, re- 
spectively, of the pulse as the most basic spatio-temporal 
pattern that sustainedly propagates in space. This prop- 
agation boundary is called dP The focus in this 
study is on the onset of excitability in active media at 
dP. 

The propagation boundary dP is determined in the pa- 



rameter space of a co-moving frame. In this frame, dP 
is caused by a saddle-node bifurcation above which trav- 
eling pulse solutions exist. As active media are usually 
based on a reaction-diffusion mechanism, this bifurcation 
is computed in a parabolic partial differential equation 
(PDE). The details of this computation can be found in 
textbooks [l7| and have been further related to the phe- 
nomenon of SD in the cortex in Ref. Q. Therefore, we 
just sketch the basic idea. The propagation boundary in 
a parabolic PDE is obtained by searching pulse profiles 
as stationary solutions in a co-moving frame. The profile 
equation in the co-moving frame must tend to the fixed 
point value of the original system for ^ — > ibcxD, with ^ 
being the spatial coordinate in the co-moving frame. A 
traveling pulse is thus equivalent to the existence of a 
homoclinic orbit satisfying a system of ordinary differen- 
tial equations, namely the profile equation system, with 
appropriate boundary conditions. 

In See. [4] we describe how the locus of the propaga- 
tion boundary dP in the parameter space is controlled 
by nonlocal spatial coupling and time-delayed feedback 
by a shift to new parameter values without adding a dis- 
tinctly new character to the spatio-temporal patterns, as 
schematically illustrated in Fig. [21 Since the propagation 
boundary is essentially a feature of an active medium 
with one spatial dimension, we can limit our main in- 
vestigation to one-dimensional systems. Yet, to get a 
picture of the patterns that arise in the neighborhood 
of dP, we will end Sec. [2] by a brief review of patterns 
in higher spatial dimensions, and provide definitions of 
weak excitability and weak susceptibility. 



2.3. Weak excitability and weak susceptibility 

The variety of qualitatively different spatio-temporal 
patterns in higher excitable regimes, i. e., beyond dP to- 
wards the oscillatory regime (to the left of dP in Fig. [T]), 
provides the foundation for a classification of spatial ex- 
citability in active media. This is in so far analogous 
to the classification into type I and II of excitable local 
elements (Sec. 12. ip as both classifications are built on 
the patterns that emerge. However, the classification of 
local excitability defines different types of mechanisms, 
whereas the classification of spatial excitability charac- 
terizes rather the degree of excitability. 

Excitability is described by ordinary differential equa- 
tions in local elements and, by partial differential equa- 
tions in spatial media. Roughly speaking, the additional 
independent spatial variable allows for more complex pat- 
terns, due to the fact that the spatial dimension renders 
the phase space infinite-dimensional. In fact, the pat- 
terns can become even more complex if spatial dimen- 
sions are increased from one to two (for example spiral 
waves in retinal spreading depression [l^) or three di- 
mensions (Winfree turbulence of scroll waves in cardiac 
fibrillation [11]). 

Close to the propagation boundary, the complexity of 
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emerging patterns is largely independent of the number 
of spatial dimensions. This fact is also paraphrased as 
weak excitability [13] ■ This term is used for active media 
to indicate that either no reentrant patterns occur (de- 
scribed by the rotor boundary dR [16[, see Figs. [T] and 
[2]) or that the rotation period is large enough, so that 
the front of the pulse does not interact with its refrac- 
tory back. The onset of interactions between front and 
refractory back in a reentrant wave pattern is described 
by the meandering boundary dM, see Fig. [TJ To the 
left of this boundary, the core of a freely rotating spiral 



dynamics of the species, called u, are coupled by diffusion 



wave performs a meandering pattern [2ll . |22| | , whereas to 



the right of dM the spiral core follows a rigidly circular 
rotation [^KFig. [21). Changing excitability parameters 
further, the spiral core can start to perform more com- 
plex manoeuvres to avoid the refractory zone (complex 
boundary dC) , sec Fig. [TJ 

Patterns of spreading depression in chicken retina are 
observed in vitro in the complex regime to the left of dC 
[l8l |. but in human brain tissue they occur close to OR. 
It was predicted that the window of cortical excitability 
lies between dR and dP [2^ and this seems to be con- 
firmed by a functional magnetic resonance imaging study 
in migraine, mapping spatio-temporal patterns of symp- 
tom reports onto the folded cortical surface [11]. These 
patterns are similar to particle- like or retracting wave 
segments that occur between dR and dP [1^, [23, [2^. 
Transient patterns can also be observed at even lower 
excitability, limited by the spreading boundary dS [29| 
(Fig. [2]). Since the regime of retracting wave segments 
is not identical with weak excitability (absence of front- 
back interactions) it was called weak susceptibility based 
on a susceptibility scale that can also be operationally 
defined in experimental systems [25l | . The control of ex- 
citability in media being weakly susceptible to pattern 
formation can be investigated in a model with one spa- 
tial dimension, because dP is essentially a property of a 
ID parabolic PDF. 



du 
'dt 



f{u) -V 
f{u) = u - —. 



d^u 
dx^ ' 



(1) 
(2) 



Here v is a bifurcation parameter. In these non- 
dimensional equations the diffusion coefficient, which 
merely scales space, is set to unity. 

Eqs. ([I|)-(lll) were proposed by Hodgkin as a model for 
spreading depression based on the bistable dynam- 
ics suggested by Grafstein [l3|. The front profile of the 
species u and its speed can be calculated analytically for 
the cubic polynomial form of ,f{u), i.,e., the generic form 
of bistable dynamics (see for example the textbook ref- 
erences [U, [12]). In the context of control one should 
think of the two stable states as one being a physiologi- 
cal (healthy) state and the other one being a pathological 
(depolarized) state, and the latter invading the former, 
which shall be prevented by control. 

As is shown in Sec. [H control of the onset of propaga- 
tion depends essentially on the interaction of the healthy 
state with the pulse front via the augmented transmission 
capability. However, in the system defined by Eqs. P])- 
^ there does not exist a boundary like dP (Sec. 12. 2p . 
In other words, there is no abrupt onset of the capacity 
to propagate fronts with finite speed. Instead, the het- 
eroclinic orbit, which corresponds to the front profile in 
a co-moving frame, exists for the whole bistable regime, 
in particular also for arbitrarily slow velocities of the co- 
moving frame including zero (standing front) [3l|. To 
obtain propagation onset behavior, we need to add in- 
hibitor dynamics, i.e., a second recovery species for the 
Hodgkin-Grafstein model. 



3.2. Extension by inhibitor dynamics 



3. REACTION-DIFFUSION WITH 
AUGMENTED TRANSMISSION CAPABILITY 

The generic framework to model spatial excitability is 
a reaction-diffusion system of activator-inhibitor type. In 
fact, already a single species model has the capacity to 
propagate fronts [30|. Such a model was originally sug- 
gested for SD by Hodgkin and Grafstein [l^] (Sec. 13. ip . 
We extend this generic model by choosing appropriate in- 
hibitor dynamics (Sec. 13. 2p to obtain pulses with the on- 
set saddle-node bifurcation dP and introduce augmented 
transmission capabilities (Sec. 13.3"]) . 



3.1. Hodgkin-Grafstein model 

A single species model of bistability has the capacity 
to propagate one state into the other if these bistable 



With a proper choice of a second species whose dy- 
namics is coupled to u, the local system can change from 
bistable to excitable dynamics, which is needed to ob- 
serve and compute dP. The natural choice is to use the 
bifurcation parameter v of the Hodgkin-Grafstein scheme 
as the additional species, namely as an inhibitor: 



dv 

m 



eg{u,v) 



(3) 



The rate equation g(w, v) of the inhibitor determines the 
dynamics in the refractory phase of the pulse, where it 
recovers the initial healthy state that was recruited into 
the pathological state of species u. In this two-species 
model, u is called the activator. Inhibitor dynamics usu- 
ally changes on a slower time scale £ ^ 1. 

There are several models in the neuroscience litera- 
ture obeying this structure in E qs. ffl-pi ), both models 
with local excitability of type I SlTssland of type II 
[s^. [35I. [36t . These models are, however, models of action 
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potentials describing the propagation of normal electro- 
physiological activity along single nerve fibers. Therefore 
these models do not relate directly to SD, because SD is 
caused by large-scale neuronal activity I37l. |38| . There- 
fore the analogy to models as in Rcfs. [311. |33|J34. Issl. [36j 
is mainly formal in the mathematical structure but not 
in its biological interpretation. For example, in the 
Hodgkin-Grafstcin Eqs. ([I])-© the activator is K'^ while 
in action potential models it is the membrane potential. 

Models exhibiting type I excitability are also capable 
to show type II behavior for a certain parameter regime 
but the reverse is not necessarily true. This suggests 
that type II behavior is more generic, because it can be 
observed in both types of models and only requires a 
rate function g{u,v) that is linear in both arguments. 
We introduce inhibitor dynamics as 

g{u,v) = u + lS — '-fv. (4) 

With this choice, the local excitable dynamics is of type 
II. The parameter /3 and 7 determine the precise con- 
figuration of trajectories in phase space that lead to ex- 
citability and are thus, in addition to e, the only two 
parameters that control excitability of type II in a nor- 
mal form model. 

Eqs. correspond to the well studied FitzHugh- 

Nagumo (FHN) systems, which is, in fact, used as a 
paradigm for excitable systems [Tsj . Pattern formation in 
a net of FitzHugh-Nagumo elements under the influence 
of time-delayed feedback was investigated with the aim 
to increase the coherence of noise-induced wave patterns 

3.3. Augmented transmission capability 

We extend the FHN model by introducing augmented 
transmission capability into Eqs. as a feedback 

loop 

( e"i9tw ) ~ ( g{u,v) ) ^ ^ ( « ) ^ ^ ( " ) ' 

where D is the local diffusion operator, and H represents 

the augmented transmission capability. In Eqs. 

we have considered diffusion in the activator species u 

only: 

The augmented transmission capability 

H = /v F (7) 

is described by the control strength K and the control 
matrix 

F = ( p"" ) , (8) 




V s 



FIG. 3: Diagrammatic view of control by augmented trans- 
mission (CAT). 



whose elements Fij are operators which represent three 
individual steps of the control by augmented transmis- 
sion (CAT) , namely (i) selecting a species j whose trans- 
mission capability is augmented, (ii) creating the control 
force from this species, and (iii) feeding this control force 
back into the dynamical variable i, see Fig. |31 

Formally, this can be represented by splitting 
Fij = AijF into the components of a coupling matrix A, 
which represents the coupling scheme, e.g., 

— — (2;). (») 

and an operator F that creates the type of control signal. 

The control signal F is generated from the input vari- 
able s = u or s = u by time-delayed or nonlocal feedback. 
For exarnple, in time-delayed feedback as introduced by 
Pyragas [4I , the operator F creates the difference between 
the time-delayed signal s{t — t) and its current counter- 
part s{t). We study this type and two nonlocal types of 
coupling. The control force F of these different types, as 
a function of the signal s, is described in detail in Sec. 14.11 
and l4.2l It should be noted that in Eq. ([5]), H represents 
cortical circuits and neurovascular coupling that is, 
internal control loops. 

The four coupling matrices in Eq. ([9]) lead to four 
transmission pathways, also termed coupling schemes: 
two self-couplings uu and vv, and two cross-couplings 
vu and uv, corresponding to the indices ij in Eq. ([8]). 
Note that other coupling schemes, for example diag onal 
coupling or coupling with a rotation matrix |42|. [4^. 
can also be investigated within this framework. 



4. SUPPRESSION OF PULSE PROPAGATION 

We perform simulations of the RD system in 
Eqs. in one spatial dimension with RD parameter 

settings (e=0.1, /3=0.85, 7=0.5) for the weakly excitable 
regime close to the propagation boundary dP (Sec. 12. 2| ). 
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The RD system is extended by three types of couphng 
in the framework of Eq. ([5]), i.e., two types of nonlocal 
spatial coupling (isotropic and anisotropic) and one type 
of time-delayed feedback. In the context of spreading 
depolarizations, nonlocal spatial and local time-delayed 
couplings represent neural structural and functional con- 
nectivity and neurovascular feedback in the cortex. This 
is discussed in more detail in Ref. For each type 
of coupling there are four principal coupling schemes, 
two self-coupling schemes (mm, vv) and two cross-coupling 
schemes {uv,vu), defined by the coupling matrix A in 
Eq. dH). 

We initialize each simulation with a stable pulse profile 
solution of the RD system in Eqs. (HJ-lH]) to investigate 
the effect of control on RD excitability. This is accom- 
plished by testing whether this specific initial condition 
of the free system (RD only) is in the basin of attrac- 
tion of the homogeneous steady state of Eq. ([5]), i. e., the 
controlled system (RD-CAT). There are two CAT param- 
eters: The control gain K in Eq. and a spatial (S) 
or temporal (r) control scale, repectively, which will be 
introduced in the following sections. For a large range of 
these two CAT parameters, we determine whether pulse 
propagation is suppressed or not. The propagation is 
suppressed if the excitation dies out, so that the system 
approaches the homogeneous steady state. In this case, 
CAT is considered successful. In the reverse case, any 
sustained spatio-temporal pattern that evolves from the 
initial conditions (free pulse solution), after the CAT is 
" switched on" , is considered an unsuccessful control since 
the activity is not completely suppressed and the homo- 
geneous steady state is not reached. 

For the following simulations we adopt an active 
medium with a spatial extension of L = 160. As spa- 
tial increment in the discretized Laplacian D we take 
Sx = 0.2. All simulations are run for 2000 time units 
and use an Euler forward algorithm with discretization 
St = 0.00125. The spatial and temporal width of the free- 
running activator pulse, Aa; and At, respectively, serve 
as reference space and time scales Q. 



in Fig. [5l For the two self-coupling schemes uu and vv, 
shown in Fig. O (a) and (d), respectively, the sign of the 
gain parameter K determines the effect of the nonlocal 
connection. Pulse propagation can only be suppressed for 
K > 0. For the two cross-coupling schemes uv and vu, 
shown in Fig. [5] (b) and (c), respectively, the sign of the 
gain parameter K changes at (5 w Ax for pulse suppres- 
sion, and these two cross-coupling schemes show similar 
control domains with respect to reflection K —K. 

The sign of K in the control domains of the self- 
couphng schemes (Fig. [5] (a) and (d)) can be qualitatively 
understood by considering the effect of nonlocal connec- 
tions upon the homogeneous steady state in the limit 
(5 — > 0. For if > 0, this limit corresponds to diffusively 
coupled elements. In general, the homogeneous steady 
state is stabilized by diffusion against small inhomoge- 
neous perturbations. In the same way, a local pertur- 
bation is leveled by a nonlocal connection in the form of 
Eq. PH)) for self-coupling. In the diffusion limit, however, 
K would increase the diffusion coefficient, which causes 
the pulses to become broader. Yet, a qualitative change 
in the dynamics cannot occur by changing the diffusion 
coefficient. Therefore the effect of suppressing pulses de- 
pends on the nonlocal character of the connection, which 
must extend at least over a distance of about 20% of the 
pulse width Ax. 

The control domains of the cross-coupling schemes 
(Fig. [5] (b)-(c)) are qualitatively different from self- 
coupling. K changes its sign &i 5 ~ Ace within one 
scheme and accordingly the situation of pulse suppres- 
sion is more complex. There is a long-range regime 
{5 > Ax) and a short-range regime {S < Ax). The sign 
of K for pulse suppression in the long-range regime de- 
pends mainly on the backward connection (solid arrow 
in Fig. [Ha)), as will be shown in the next section. In 
the short-range regime, the effect of the forward connec- 
tion (dotted arrow in Fig. HJa)) seems to dominate and 
reverse the effect observed in the long-range regime. 



4.1. Pulse suppression by nonlocal spatial coupling 

In this section we present results on two types of spatial 
coupling. 



4.1.1. Isotropic coupling 




space X 



Nonlocal isotropic spatial coupling is defined as 

F{s) = s{x + 5,t) + s{x -S,t)- 2s{x, t) (10) 

where s = u ox s = v. 

Regimes in which reaction-diffusion pulses are sup- 
pressed by additional nonlocal connections are calculated 
for the parameters 5 and K and represented by gray areas 



FIG. 4: (Color online) Schematic diagram illustrating nonlo- 
cal spatial coupling and time-delayed feedback, (a) Nonlocal 
spatial coupling: isotropic coupling (dashed and solid arrow), 
and anisotropic backward coupling (solid arrow), (b) Time- 
delayed feedback: backward coupling of traveling pulses. 
Pulse profiles are shown in red (solid) and green (dashed) 
for the activator and inhibitor, respectively. Dotted: time- 
delayed profiles. 
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FIG. 5: Control planes of isotropic spatial coupling spanned 
by the two CAT parameters: control gain factor K and non- 
local space scale 5 normalized to pulse width Aa; (left scale) 
and in spatial units (right scale), (a) Activator self-coupling 
scheme uu, (b) cross- couplings uv (inhibitor signal fed back to 
activator rate equation) and (c) vu (reverse) , and (d) inhibitor 
self-coupling vv. Suppression of pulse propagation is marked 
by gray control domains. Parameters: e — 0.1, /3 = 0.85, 
7 = 0.5. Aa; = 8.65. 



4-1.2. Anisotropic backward coupling 

Nonlocal connections can also be introduced in only 
one direction 



F{s) =s(a;±(5,t) -s(a;,t). 



(11) 



This directed connectivity would correspond to 
anisotropic nonlocal coupling in a two dimensional 
excitable medium. One reason to investigate this type 
of connectivity is to obtain a better understanding of 
the results in the previous section by separating effects 
of forward and backward connections (Fig. [4] (a)). 
Moreover, the functional and structural connectivity of 
the cortex is realistically modeled as an anisotropic (and 
also inhomogeneous) medium due to the patchy nature 
of nonlocal horizontal cortical connections. While 
anisotropics in the cortical connections are usually 
considered to merely cause variations in wave speed in 
different directions, inhomogeneities are known to cause 
wave propagation failure [44| . In contrast, our focus 
is on the change in excitability by anisotropic nonlocal 
coupling that leads to suppression of wave propagation. 

We limit your investigation to the backward connec- 
tion. This corresponds to the plus (minus) sign in 
Eq. (jlip for pulses propagating in the positive (nega- 
tive) X direction (solid arrow in Fig. Ufa)). We choose 
the backward connection because this type of coupling 
shares many properties with local time-delayed feedback 
coupling, as we will show by comparing the results from 
this section with those of the next section (Sec. 



The control domain of the uu self-coupling scheme 
(Fig. [S] (a)) is very similar to the isotropic one shown 
in Fig. [5] (a). This indicates that in the uu scheme 
of isotropic coupling the backward connection accounts 
for the main contribution to suppression of wave prop- 
agation. However, we want to note that both control 
types (isotropic and anisotropic) with the same coupling 
scheme uu can differ in their efficiency within some parts 
of the gray control domain (not shown). The efficiency 
refers to the length a pulse travels after the connectivity 
is "switched on". This transient effect was investigated 
in detail in 0, [2^ , whereas in this study we concentrate 
on an understanding of the size and shape of control do- 
mains for the different schemes and types of coupling 
(Figs.EHll). 

Pulse propagation is not suppressed by the vv self- 
coupling scheme for anisotropic backward coupling for 
any control parameter within the investigated range 
(Fig. [H^d)). This, however, does not mean that the pulse 
propagates essentially unchanged. Any sustained spatio- 
temporal pattern other than the homogeneous steady 
state that evolves from the free initial pulse solution — 
after the connectivity is "switched on" — is considered as 
unsuccessful pulse suppression. In fact, for the vv scheme 
of anisotropic coupling in the expected control domain of 
corresponding Fig. O (d), the homogeneous steady state 
becomes unstable and stationary spatial patterns emerge 
after the initial pulse is suppressed. 

The control domains of the cross-coupling schemes 
(Fig. [6] (b)-(c)) arc similar to Fig. [5] (b)-(c) in the long- 
range regime, except for a failure of suppression in the 
isotropic uv scheme for K > 0.2. The long-range regime 
for anisotropic backward cross-coupling extends over the 
nearly complete spatial scale of i5 (except for very small 
5). Thus, K does not change its sign. Also, there is a 
perfect reflection symmetric pattern in the control planes 
with respect to the axis i^T = for the two cross-coupling 



3 
2.5 

2 
1.5 

1 

0.5 


3 
2.5 

2 
1.5 

1 

0.5 




(a) 



I 



(b) 


uv 







(c) 



(d) 



-1 -0.5 



0.5 1 -1 -0.5 



0.5 1 



25 

20 

15 

10 

5 



25 

20 

15 

10 

5 





FIG. 6: Control planes of anisotropic backward spatial cou- 
pling. Parameters and notation as in Fig. [S] 
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FIG. 7; Control planes of time-delayed feedback coupling 
spanned by the two CAT parameters: control gain factor K 
and time delay r normalized to pulse width At (left scale) 
and in temporal units (right scale), (a) Activator self-coupling 
scheme uu, (b) cross- couplings uv (inhibitor signal fed back to 
activator rate equation) and (c) vu (reverse) , and (d) inhibitor 
self-coupling vv. Suppression of pulse propagation is marked 
by gray control domains. The black lines denote saddle-node 
bifurcations. The dashed lines mark the values of control 
parameters for which the u-amplitude of a single local FHN 
system, stimulated in u by I{t) (shown as inset in (d)), is 
decreased by 10% compared to the u-amplitude of the uncon- 
trolled system. Parameters as in Fig. [5] At — 10.73. 



schemes. 

The reason for the reflection symmetry in the pattern 
of control domains for cross-coupling schemes and, fur- 
thermore, the explanation for the observed signs of the 
gain parameter K for the control domains of all coupling 
schemes for anisotropic backward coupling arc given in 
Sec. 14.31 together with the explanation of some of the re- 
sults we have obtained for local time-delayed feedback 
control (Pyragas feedback), which will be described in 
the next section. 



4.2. Pulse suppression by local time-delayed 
feedback 



In case of local time-delayed feedback the control force 
F is given by 



F{s) = s{x, t — t) ~ s{x, t), 



(12) 



where r is the delay time. Note the formal similarity to 
Eq. (fTT|) . This control method was first introduced by 
Pyragas for chaos control. In Fig. [3](b), this control 
method is illustrated for the uu coupling scheme: at each 
spatial location x the activator u at time t receives the 
signal from the same location but at the past time t — t. 
That means particularly for the dynamics of the front 



that the deviation from the homogeneous fixed point is 
fed back. 

The domains of successful control, i.e., pulse suppres- 
sion, in the (A', r) plane are shown as gray areas in Fig.[7l 
For successful control, the pulse dies out and the system 
returns to the homogeneous steady state, as shown ex- 
emplarily in the space-time plot of Fig. [8] (b). 

For time-delayed feedback a domain of successful con- 
trol exists for each coupling scheme. As for the spatial 
coupling schemes in Fig. [5] and [H the two cross-coupling 
schemes uv and vu show a reflection symmetry with re- 
spect to the axis K — Q. This symmetry is explained 
by the effective parameters introduced in Sec. 14.31 Cou- 
pling schemes that feed back the signal to the activator 
can suppress pulses for positive K . For coupling schemes 
that feed back the signal to the inhibitor successful con- 
trol is possible for negative values of K. This is simi- 
lar to the results of anisotropic backward coupling (cf. 
Fig. [5]). This similarity is due to the similarity of the 
signals the pulse receives through the feedback. In both 
cases the wave front gets the feedback from the homo- 
geneous steady state which is ahead of the wave, in a 
temporal or spatial sense, respectively. 

The main differences between the control domains of 
time-delayed feedback and spatial anisotropic coupling 
are that in case of time-delayed feedback the control do- 
mains are limited towards large values of t (upper bor- 
ders of gray domains in Fig.[7](a-d)), and, moreover, that 
there is a control domain also for wd, the inhibitor self- 
coupling scheme. The limitation of the control domains 
towards large values of r is caused by the formation of 
tracking patterns [i^ (Fig. [5]), which is further explained 
in Sec. 14.2.21 These patterns emerge by delay- induced os- 
cillations. The local dynamics becomes bistable due to a 
saddle-node bifurcation (black solid lines. Fig. [7|). This 
is demonstrated in the next section. 



4-2.1. Effect of time-delayed feedback in a local excitable 
system 

In this section we investigate the effect of time-delayed 
feedback on a single local excitable element. There are 
two main effects. First, depending on the sign of the con- 
trol parameter K, the activator amplitude is reduced, 
which selects the location of the control domain in the 
control plane with respect to the axis K = (Fig. [7]). 
Furthermore, we show that a minimum amplitude reduc- 
tion of about 10% is needed to suppress pulse propaga- 
tion when these elements form an active medium (lower 
bounds of the control domains). Second, for too large a 
value of T, the local dynamics becomes bistable (upper 
bounds of the control domains). 

To obtain qualitative insight into how time-delayed 
feedback operates, we perform numerical simulations in 
order to compare trajectories starting from the same ini- 
tial conditions with and without time-delayed feedback. 
Figure O shows exemplary superthrcshold phase space ex- 
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FIG. 8: Pattern formation for e = 0.1, /3 = 0.85 and 7 = 0.5: 
(a): pulse propagation without control; (b): suppressed pulse 
via nonlocal coupling (K = 0.2, 5 = O.58A2:). The horizontal 
line at t = 25 marks the onset of control. The distance that 
the transient excitation spreads after control is activated be- 
fore the medium relaxes to the homogeneous fixed defines the 
transient activation radius (TAR) or tissue at risk, (c): delay- 
induced oscillatory tracking pattern {K — 0.5, r = 2.8At, 
onset of control at t — 21) 



cursions with and without (A' = 0) time-delayed feed- 
back (solid and dashed trajectories, respectively). The 
system is initialized in the interval [tQ—T;tQ[ (history 
function) with the fixed point value {u — u* and v = v*) 
and for with a superthreshold value {u = u* + 0.5 and 

V = V*). 

For coupling schemes that feed back into the activator 
rate equation, i. e., uu and uv, the control force acts par- 
allel to the u axis. In Fig.[n](a) the trajectories with and 
without control for the uu coupling scheme are plotted. 
The direction of the control force is denoted by the hori- 



zontal arrow. In order to reduce the amplitude in u, the 
control force has to be directed towards the fixed point, 
i.e. it has to be negative for u > u* and v > v* . This is 
the case for t < r if > 0, because the history function 
is initialized as {u*,v*). 

For the vu and vv coupling scheme the control force 
acts parallel to the w-axis. In Fig. [9] (b) the trajectories 
with and without control are shown for vv coupling. The 
direction of the control force is denoted by the vertical 
arrow. To get lower amplitudes of u the control force has 
to be directed towards the opposite direction of the fixed 
point. This is the case for t < r if K < 0. 

In the following the reduction of amplitude u of a sin- 
gle local FHN system is investigated in order to obtain 
the lower boundaries (dashed) of the control domains in 
Fig. [71 To excite the system a stimulation current I{t) 
(inset in Fig. [7]^d)) is added to the activator u. Choos- 
ing for each coupling scheme the proper sign of control 
gain K to reduce the activator amplitude, time-delayed 
feedback is activated and the maximum amplitude Umax 
is observed. By performing a bisection method for each 
value of K, the proper r is detected that reduces Umax by 
10% from the original amplitude after stimulation with- 
out feedback {K = 0). The dashed lines in Fig. [7] mark 
the positions where the it-amplitude is reduced by 10% 
for a fixed stimulation. These lines form the lower bound- 
aries of the control domains. 

Also the upper boundaries of the control domains can 
be understood by investigating the single system with 
feedback. For large K and r the system becomes bistable 
in the sense that in addition to the stable fixed point a 
stable limit cycle occurs. This limit cycle appears in a 
saddle-node bifurcation of limit cycles, as was also found 
in a system of two delay-coupled FHN-systems [i^ \^ . 
The limit cycle emerges if (i) t is sufficient large to allow 
for recovery to the fixed point, and (ii) the feedback and 
hence the gain parameter K is strong enough to push the 
system beyond the threshold. 

Again with the help of a bisection method the saddle- 
node bifurcation lines for the single system are deter- 
mined in the K-t plane. In Fig. [7] they are plotted as 
black solid lines. These lines mark the boundaries where 
the local dynamics becomes bistable and control fails. 
Only in case of uu coupling the boundary of the con- 
trol domain deviates appreciably from this bifurcation 
line. This is due to the stabilizing effect of diffusion that 
damps out local oscillations. Since in the interspace the 
single system is bistable and thus the medium is able 
to perform homogeneous oscillations, the difference be- 
tween the upper boundary of the control domain and the 
bifurcation line depends on the initial conditions that are 
chosen to locally stimulate the medium. 

Investigating a single FHN-element with time-delayed 
feedback provides a qualitative understanding of the dy- 
namics of the controlled spatial system: The sign of K 
and the form of the control domains can be understood. 
Thus the control of local excitability provides control of 
the global spatial excitability. This will be quantified 
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in Sec. 14.31 In the next section, we present the spatio- 
temporal patterns that emerge above the upper bound- 
ary of the control domain, i.e., the bistable domain of the 
single local element. 

4-2.2. Delay-induced oscillatory pattern formation 

As we have seen in the previous section, the local FHN 
dynamics with time-delayed feedback becomes bistable 
for too large values of K and r. Due to this local bista- 
bility, in the spatially extended system delay induced os- 
cillatory pattern formation occurs (Fig. (I5]c)): After a 
local stimulation the medium performs local oscillations 
that spread slowly: Within each local oscillation the exci- 
tation spreads a short way, and thus gradually new parts 
in the neighborhood become excited. The area of oscil- 
latory excitation grows slowly all over the medium. We 
have observed these patterns for each coupling scheme 
beyond the upper boundary of the control domain. 

This effect will again play an important role in sec- 
tion |43?2l that focuses on the propagation boundaries in 
(e,/?) - parameter space. 




FIG. 9: Phase portraits of a local FHN system with (solid tra- 
jectory) and without (dashed trajectory) time-delayed feed- 
back control, (a) uu coupling: K=0.03 (b) vv coupling: K=- 
0.7. Parameters: e — 0.1, (3 = 0.85, 7 — 0.5, r = 5, no 
diffusion. The solid arrow denotes the action of the control 
force. 



4.3. Description through effective parameters 

In this section we further invastigatc the two control 
schemes of spatial anisotropic backward coupling and lo- 
cal time-delayed feedback to achieve an analytic approx- 
imation. In Sec. 14.3.1] we use this to describe the change 
of excitability and in Sec. 14.3.21 we investigate how this 
is reflected in the parameter space. 

For our approximation we assume that the front dy- 
namics governs the location of the propagation boundary 
dP. Therefore we focus on the part of control that acts 
on the front dynamics. For the coupling types of time- 
delayed feedback and nonlocal anisotropic coupling the 
front of the pulse receives the feedback of the homoge- 
neous steady state. The main idea is to replace the time- 
delayed (in Eq. ^) or the space shifted (in Eq. (|lip ) 
quantities by their fixed point values {u{x, t ^ t) ~ 
u(x + 5,1) = u* and v{x, t ~ t) = v{x -\- 6,t) = v*). This 
transforms the controlled system, after some simple al- 
gebraic manipulations, into the form of the uncontrolled 
system (K=0), but with new effective values of the pa- 
rameters £{K), f3{K) and j{K). 

For each coupling scheme the obtained effective param- 
eters and the corresponding transformations arc shown in 
the Table |T1 Since the delayed or shifted quantities are 
replaced by the fixed point values, the obtained effective 
parameters e, (3 and 7 only depend on K and not on t 
or S. 

Although the transformations look rather complex 
there is a common feature for vv and the two cross- 
coupling schemes vu and uv: The change of the parame- 
ters is such that the fixed point value of the system does 
not change, i.e., for these three couphng types the fixed 
point after the transformation to the system with effec- 
tive parameters has the same homogeneous fixed point 
as the system without control. Since for the uu coupling 
scheme the transformation is more complex, the fixed 
point in the transformed equations differs from that in 
the original equations, i.e., the intersection of the nuU- 
clines changes its relative position with respect to the cu- 
bic nullcline. Note that this does not mean that the fixed 
point of the original system changes by adding control. 
Due to the non-invasivity of all investigated control types 
the fixed point does not change. However, the transfor- 
mation that converts the original system with control into 
a new system without control but effective parameters, 
is in case of uu coupling not invariant for the fixed point 
iu*,v*). 

The effects of parameter changes in the FHN system 
are well known 0, [2^ . As rule of thumb one can keep in 
mind that the larger e,(3 or 7, the lower the excitability 
of the system. 

In the following the K dependence is discussed: 
In case of vv coupling e does not change. For K < only 
(3 becomes larger while 7 decreases. However, the change 
of /3 dominates, i.e., the excitability decreases for K < 0. 

For vu coupling e and /3 increase for isT < 0, while 
7 decreases. In that case the influence of changing e 
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dominates, i.e., excitability decreases for if < 0. 

For uu and uv coupling e and (3 increase for if > 0, 
while 7 decreases. The influence of e is dominating, and 
therefore the excitability decreases for positive values of 
if. 

For the vv and vu schemes the inhibitor receives a feed- 
back. In these cases by simply rearranging the inhibitor 
equation the original form without control, but with the 
new effective parameters, can be obtained, and therefore 
no transformation in time and space or of u and v vari- 
able has to be performed. For the two coupling schemes, 
for which the activator receives a feedback [uu and uv) 
the equations need to be transformed in order to retain 
the original form without additional control force. In 
case of uv only the v variable is transformed, whereas 
in case of uu coupling both dynamic variables u and v 
and also time and space are transformed. However, this 
does not change the qualitative dynamics of the system, 
since transformations in u and v variable and in time and 
space correspond to rescaling only. 

For the two cross-coupling schemes the effective pa- 
rameters are symmetric with respect to if — ^ —if. This 
symmetry of the cross-coupling schemes was observed in 
the control domains of all control types (Figs. O El and 

[21). 

4-3.1. Change of excitability 

In this section we clarify the influence of the control 
parameter if on the excitability of the system, within the 
approximation of effective parameters. In particular the 
influence of [3 and 7 is investigated, since they change 
for all control schemes in opposite ways. For vv coupling 
the influence of f3 dominates and for the other coupling 
schemes the influence of e is decisive. 

For the three coupling schemes vv, vu, and uv the 
transformations from the original system to that with ef- 
fective parameters leave the fixed point invariant, i.e. the 
fixed point values with respect to the transformed vari- 
ables u and V are the same as with respect to the original 



one {u and v). Since the fixed point depends only on 
f3 and 7, the condition that the fixed point remains the 
same forces the effective parameters /3 and 7 to change si- 
multaneously. Thus they are not independent from each 
other. The parameters — /3 and i define the intersection 
with the u axis and the slope of the ii-nuUclinc of the 
FHN system, respectively. Thus, for a smaller slope of 
the w-nuUcline and unchanged fixed point, 7 has to be 
increased and (3 decreased and vice versa. For the three 
cases of uv, vu and vv the dependence between /? and 7 
yields the condition 

/3(7) = /3+(7-7)«*- (13) 

It is not intuitively clear how excitability changes. In- 
creasing 7 while respecting the invariant fixed point con- 
dition given by Ea. ([T^ . decreases f3, which yields two 
opposing effects on excitability. In the following it is 
shown that under the condition in Eq. (I13p the infiuence 
of f3 dominates. 

To investigate the influence of simultaneously changing 
f3 and 7 with unchanged fixed point, numerical simula- 
tions of a single FHN systems with effective parameters 
are performed under the condition in Eq. (|13p . To excite 
the system the stimulation current I{t) (Inset in FigfTOl) 
is added to the activator u. The response of the ac- 
tivator u is determined in dependence of 7 (if). Note 
that these simulations are done without feedback. In 
order to deduce the infiuence of f3 and 7 upon excitabil- 
ity, e = e = 0.1 is not changed during the simulations. 
Thus the chosen parameters are equivalent to the effec- 
tive parameter set of vv coupling. For the other system 
parameters (3 = 0.85 and 7 = 0.5 are chosen. 

In Fig. (TUl the maximum amplitude Umax of the acti- 
vator is plotted versus 7(if). For 7(if) < —0.05 the re- 
sponse of the system to the stimulation I(t) has a small 
amplitude. For 7(if ) > —0.05 the amplitude of the sys- 
tem blows up in a very tiny parameter range. After this 
blow-up the amplitude remains for all 7(if ) > —0.05 at 
about the same level. 

This blow-up resembles the canard explosion of the 
FHN system near the Hopf bifurcation. The difference 
is that a canard transition usually characterizes the fast 
blow-up of limit cycles, whereas in case of an excitable 
system a limit cycle does not exist, but only its ghost is 
visible from the excited trajectories. However, the result 
is that increasing 7(if ) with simultaneously adjusting P 
according to Eq. (|13p increases the response to a stimu- 
lus. We conclude that excitability increases. Therefore 
it follows that the infiuence of decreasing (3 dominates. 

Considering the ranges of 7(if ) that result for the dif- 
ferent coupling schemes by varying if in the interval 
[—1,1] the impact of each effective parameter on the cou- 
pling schemes can be determined. 

For the vv coupling scheme the 7(if ) is in the interval 
[-0.5,1.5] for if e [-1,1]. The excitation of the loop 
leads to the canard like blow-up close to 7(if ) = —0.05, 
which is equivalent to Kyy = —0.55. Hence we estimate 
for Kyy = —0.55 a change of excitability. This is in 
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FIG. 10: Maximum amplitude Umax of activator u after 
stimulation with I(t) (inset: solid line). System parameters: 
e=0.1, I3=P{K), ^—^{K). The bars below the a;-axis mark 
the ranges for the effective parameter 'y{K) for vv coupling 
and the two cross-coupling schemes uv and vu. K is varied 
in the range [—1, 1]. 

good agreement with the onset of successful control for 
time-delayed feedback in the spatially extended system 
(cf. Fig. Hid)). 

For the two cross-coupling schemes uv and vu the 
change of e dominates. jiK) lies in the interval [0.25, oo[. 
The values of 7 for these two schemes are beyond the 
canard-like transition, in the region of large amplitudes. 
For the investigated values of K the system does not 
undergo the canard-like transition from large to small 
amplitudes and hence the influence of 7 and (3 on the 
excitability is small. Therefore the influence of e is de- 
cisive. To verify this assumption we estimate the value 
of K^^/^u which is necessary to reach the propagation 
boundary dP by only considering e. The propagation 
boundary dP for (3 = 0.85 and 7 = 0.5 is reached for 
£ = 0.1123. For larger e a stable pulse does not ex- 
ists. Since the simulations were performed for e = 0.1, 
one can compute the values of K^v/vu needed to move 
from £ = 0.1 to £ ^ 0.1123, neglecting the influence of 
(3 and 7. In the to cross-coupling schemes this yields 
Kuv/vu = ±0.123, which is in very good agreement with 
the results from the simulations, where successful sup- 
pression is found for |J!r„i,/^„| > 0.12. 

Performing the same calculation for the uu coupling 
assuming that also in that case the influence of e dom- 
inates, £ = 0.1123 is reached for K = 0.056, whereas in 
the simulations already for K > 0.03 successful control 
was observed, which is still of the correct order. 

4-3.2. Shift of propagation boundary 

Above we have introduced effective parameters by re- 
placing the time-delayed or space-shifted quantities in 
Eq. and Eq. pT|) respectively, and have investigated 
the influence on excitability for one set of system param- 
eters (e = 0.1, P = 0.85, 7 = 0.5) and variable control 
parameter (K, S and r). In this section the influence 
of the effective parameters on the propagation boundary 



dP in (e, /?)-parametcr space (for 7 = 0.5) is compared to 
the one obtained through simulations with time-delayed 
feedback. Therefore \K\ = 0.2 and t = 0.5At is chosen, 
where At is the temporal pulse width of the activator u. 
The sign of the control parameter K is chosen such that 
the excitability is decreased, namely negative for vv and 
vu and positive for uu and uv. 

The propagation boundaries are obtained in two differ- 
ent ways. Those with time-delayed feedback are obtained 
by simulating the full equations and performing a bisec- 
tion method. Those propagation boundaries with effec- 
tive parameters are obtained by continuation of homo- 
clinic orbits (stationary pulse profiles) in the co-moving 
frame 0,|4|. 

In Fig. [11] the propagation boundaries dP are shown. 
Those obtained with effective parameters are the full 
lines separating different colors. The white dashed curves 
display the propagation boundaries for the full system 
with time-delayed feedback. The boundaries are in good 
agreement for vv coupling and the cross-coupling schemes 
uv and vu. For uu coupling the two propagation bound- 
ary differs from each other. The reason is that the effec- 
tive parameters are only exact for large r. 

The inset shows propagation boundaries of the uu cou- 
pling. The full line represents that obtained by effective 
parameters. The different thin lines represent the prop- 
agation boundaries obtained by simulating the full sys- 
tem with different time delays. In the range of small 
£ for increasing r the displayed propagation boundaries 
of uu coupling converge towards the propagation bound- 
ary of the effective parameters. The propagation bound- 
ary with r = 30At perfectly fits the one obtained with 
effective parameters for £ < 0.04. However, for larger 
£, this boundary (dotted) diverges sharply from the one 




£ 

FIG. 11: Shift of propagation boundary dP: full lines: ef- 
fective parameters, dashed lines: with time-delayed feedback 
r — 0.5At. Inset: propagation boundary for uu coupling, full 
line: with effective parameters, dashed lines: with different 
delay times r — 5At, T=10At, r=15At, T—30At (r is given 
in units of the pulse width At). Parameters: 7 = 0.5, K = 0.2 
for uu and uv, K = —0.2 for vu and vv. 
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with effective parameters. Tfiis is eaused by tlie delay- 
induced bifurcation that occurs for large K and r, see 
Section 14.2.11 Fig. [8] (c) shows the originating patterns 
that are already described in Section 14.2.21 After a lo- 
cal stimulation the time-delayed feedback suppresses the 
emerging pulse, as predicted by the effective parameters. 
But another effect occurs: the feedback is strong enough 
to create a delay-induced oscillation. This is due to the 
bistability of the local system, where in addition to the 
stable fixed point a stable and an unstable limit cycle 
are born in a saddle-node bifurcation. Each of the exci- 
tations can spread for a small distance until it is again 
suppressed by the feedback. Step by step this pattern 
propagates and grows. 



5. DISCUSSION AND CONCLUSION 

The increasing interest in both computational investi- 
gations of brain functioning [i^ and control of complex 
dynamics [l| has led us to combine methodologies and 
concepts from both fields to investigate the pathogenesis 
and potential treatment of brain disorders [50| . This is- 
sue sets the context for our studies. However, we expect 
that our results on controlling traveling pulses can find 
also applications in other fields of biomedical engineering, 
since traveling pulses occur in many biological systems 
[Sll. [5^. Furthermore, in this study we consider models 
and methods of generic type, i.e., on the one side, we 
study pattern formation in reaction-diffusion systems of 
activator-inhibitor type, namely the FitzHugh-Nagumo 
system. On the other side stands a universal method 
of chaos control, that is, time-delayed feedback (Pyragas 
control ) , which is used to control the reaction-diffusion 
patterns. In addition, we consider nonlocal spatial cou- 
pling as a control method. Nonlocal spatial coupling and 
time-delayed feedback have, as we have shown in Sec. IH 
a common mechanism that underlies the control of trav- 
eling pulses. 

The discussion will focus on two issues. Firstly, the 
comparison of time-delayed feedback and nonlocal spatial 
coupling, in particular the predictive power of the effec- 
tive parameters that can be introduced in the same man- 
ner in both cases (Sec. 15. 1[) . Secondly, the time-delayed 
feedback and nonlocal spatial coupling will be consid- 
ered from a broader perspective as augmented transmis- 
sion capabilities in reaction-diffusion systems with par- 
ticular respect to the corresponding cortical structures 
(Sec. (EH). 



5.1. Time-delayed feedback and nonlocal spatial 
coupling 

The control force F{s) in Eq. pT|) . i. e., the force in the 
anisotropic nonlocal type of coupling with backward con- 
nections, can be directly compared to the control force 
F{s) in Eq. (fT^ . which is time-delayed feedback (Pyra- 



gas control) applied to each element in the active media 
locally. By going from Eq. (fT^ to Eq. (|lip nonlocal con- 
nections are introduced simply by changing the position 
of the shift operator from the first to the second argu- 
ment of signal s{x,t). In other words, Pyragas control is 
translated from the temporal to the spatial domain, as 
illustrated in Fig.Ul This is compatible because the pulse 
is stationary in the co-moving frame and thus the speed 
of the pulse relates space to time scales. If 5 and r are 
normalized to the pulse width in the spatial and tempo- 
ral domain. Ax and At, respectively, the common effect 
of nonlocal coupling and time-delayed feedback on trav- 
eling pulses is reflected in similar locations of the control 
domains in the control planes in Fig. [5] and [T] 

The analogy between nonlocal coupling and time- 
delayed feedback, of course, oversimplifies the situation. 
However, the use of both types of coupling has been pro- 
posed for control of spatio-temporal chaos in spatially 
extended systems based on the idea of stabilization of 
unsta ble p eriodic patterns embedded in spatio-temporal 
chaos |53| . Using both types of coupling simultaneously, 
it has been demonstrated through numerical analysis 
that unstable roll patterns in a transversely extended 
three-level laser model can be stabilized. The motivation 
for this combined approach to control unstable periodic 
patterns lies in the noninvasive character of both non- 
local coupling and time-delayed feedback if the unstable 
periodic patterns is approached. 

Our motivation to study nonlocal coupling and time- 
delayed feedback is somewhat different from that of chaos 
control, because we do not want to stabilize unstable 
periodic patterns. We investigate the control of trav- 
eling pulses by comparing nonlocal coupling with time- 
delayed feedback using both types separately. We sug- 
gest a method to predict the effect of the gain factor 
K on excitability by identifying the controlled system 
with the free system with effective parameters, based on 
the idea that the effect of control makes its main con- 
tribution upon the front dynamics. Under this assump- 
tion, we suggest to replace in both Eq. pT|) and Eq. 
the shifted quantities by their fixed point values, i.e., 
s(a;, t — t) ~ s{x + 6,t) = s* , where the signal s was cho- 
sen to be either u or v and s* is the corresponding fixed 
point value. 

Despite its simplicity, this method yields accurate pre- 
dictions for traveling pulses invading the homogeneous 
steady state. The resulting equations of effective param- 
eters, given in Tab. [H describe by simple algebraic rela- 
tions how control changes excitability (Fig. [2]). However, 
in which direction excitability changes with changing K 
in the four coupling schemes has to be investigated by 
evaluating the combined effect on all three effective pa- 
rameters e, P and 7 (Fig. [TU|) . 

The noteworthy feature of our method is that the ef- 
fective parameters describe a shift in the excitability of 
the local elements. By investigating time-delayed feed- 
back in a single local element it is shown that the onset of 
pulse suppression as well as the boundedness of the con- 
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trol domains for large K and r can be explained by the 
local dynamics. Furthermore, transferring the results of 
time-delayed feedback to the coupling type of anisotropic 
nonlocal backward coupling one finds the same depen- 
dences. 

5.2. Reaction-diffusion with augmented 
transmission as a hybrid model for the cortex 

Control introduces augmented transmission capabili- 
ties in the reaction-diffusion model (Fig. [3]). The result- 
ing hybrid model in Eq. ([5]) combines the two major sig- 
nalling systems in the brain, namely local coupling by 
diffusion, termed volume transmission, and nonlocal cou- 
pling in the spatial domain described by Eqs. (|10p and 
pT|) or in the temporal domain by Eq. p2)) . The aug- 
mented transmission capabilities are typical for synaptic 
transmission and neurovascular coupling. For example, 
the change of sign in the gain parameter K for pulse 
suppression (Fig. [5] b,c) is reminiscent of the Mexican- 
hat type functional and structural connectivity pattern 
in the cortex. Also time-delays of the order of sec- 
onds, that is, the order of the width of the pulse profile 
At in spreading depolarizations, can occur in synaptic 
transmission if metabotropic ion channels are involved, 
like metabotropic glutamate receptors, which have in- 
creased open probabilities in the range of seconds after 
their activation. Moreover, typical latencies of this or- 
der result from the neurovascular coupling. Therefore, 
the augmented transmission capabilities represent inter- 
nal neural circuity that is complementary to the vol- 
ume transmission introduced in the original Hodgkin- 
Grafstein Eqs. (H])-© as a model for spreading depolar- 
izations. 

The first hybrid model for spreading depolarizations 
that also combines the two major signalling systems in 
the brain has been studied by Reggia and Montgomery 
[53 |. Potassium dynamics was modeled by a quadratic 
rate function (cf. Eq. ([2])) and coupled to a neural net- 



work that mimics cortical dynamics and sensory map 
organization. The key result of this model is that at 
the leading edge of the potassium pulse, the elsewhere 
largely uniform neural activity was replaced by a pat- 
tern of small, irregular patches and lines of highly active 
elements, which could explain neurological symptoms in 
migraine. However, there is no feedback of the network 
activity to the potassium reaction-diffusion pulse. There- 
fore this hybrid models can not address the questions of 
the controversially discussed Hodgkin-Grafstein mecha- 
nism [5^, [5^ because the augmented transmission ca- 
pabilities work only one way. In a more general con- 
text, such hybrid models are similar to creation of spatio- 
temporal networks in addressable excitable media, which 
are studied in the chemical Bclousov-Zhabotinskii reac- 
tion [57j . 

Volume transmission described by the Hodgkin- 
Grafstein mechanism was long thought to be the main 
factor that causes the propagation in spreading depo- 
larizations [s^ although synaptic transmission and gap 
junction coupling were suggested to provide an alterna- 
tive mechanism (56. [s^. Hybrid models can address these 
controversial issue and may help to provide insight into 
the spread and control of pathological pulses in the brain. 
Our emphasis is on understanding the role of internal cor- 
tical circuits that provide augmented transmission capa- 
bilities and can prevent such events. However, a long- 
term biomedical engineering therapeutic aim is also to 
design strategics that either support the internal cortical 
control or mimic its behavior by external control loops 
and translate these methods into applications. 
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